Create function to check rasters against ArcMap outputs

find_diffs <- function(test_rast, key_rast) {
  ### normalize both rasters to 1s. 
  ### Presumably, extents, CRS, and res should be identical, so 
  ### stack should work just fine:
  s_rast_norm <- test_rast / test_rast
  k_rast_norm <- key_rast / key_rast
  tmp_stack <- stack(s_rast_norm, -k_rast_norm)

  diff_rast <- calc(tmp_stack, fun = sum, na.rm = TRUE)
  diff_rast[values(diff_rast) == 0] <- NA
  return(diff_rast)
}

plot_diffs <- function(test_rast, key_rast, title = NULL, pal = 'Green-Orange') {
  tmp <- find_diffs(test_rast, key_rast)
  if(is.null(title)) title <- names(test_rast)
  plot(tmp, main = title, col = hcl.colors(3, pal))
}

check_rast <- function(test_rast, key_rast, tol = .002) {
  n_cells_key <- sum(!is.na(values(key_rast)))
  
  cellcount <- sum(!is.na(values(test_rast)) - n_cells_key)
  cellcount_check <- cellcount / n_cells_key %>%
    abs()
  
  ### find cells where the two rasters disagree 
  diff_rast <- find_diffs(test_rast, key_rast) %>%
    abs()
  
  celldiff <- sum(values(diff_rast), na.rm = TRUE)
  celldiff_check <- celldiff / n_cells_key
  
  check <- celldiff_check <= tol & cellcount_check <= tol
  
  return(data.frame(good      = check,
                    diff_tot  = celldiff,
                    celldiff  = celldiff_check, 
                    cellcount = cellcount_check))
}

Create wind mask

Here we can simply use the raster::reclassify() function just as we did in ArcMap model builder. We’ll need to set up a matrix of original values and the new values. There’s also another way to do the same thing, shown after reclassify().

wind_raw_rast <- raster('data/Wind2002.tif')
plot(wind_raw_rast, main = 'Raw wind raster', legend = FALSE)

reclass_mtx <- matrix(data = c(c( 1,  2, 3, 4, 5, 6, 7),
                               c(NA, NA, 1, 1, 1, 1, 1)),
                      ncol = 2, byrow = FALSE)
reclass_mtx

##      [,1] [,2]
## [1,]    1   NA
## [2,]    2   NA
## [3,]    3    1
## [4,]    4    1
## [5,]    5    1
## [6,]    6    1
## [7,]    7    1

wind_mask_rast <- raster::reclassify(wind_raw_rast, rcl = reclass_mtx)
plot(wind_mask_rast, main = 'Reclassified wind raster', legend = FALSE)

But we can also use vector indexing to directly replace the values too!

### first make a copy of the wind_raw_rast so we can change it
wind_mask2 <- wind_raw_rast

### turn low wind values to NA:
values(wind_mask2)[values(wind_mask2) < 3] <- NA

### note, we haven't changed the other values to 1 yet; technically
### that's fine - a mask really only needs to have data cells and
### NA cells.  But let's change them anyway - here's a trick I like to
### use to quickly set all non-NAs (and non-zeros) to 1: simply 
### divide the raster by itself!
wind_mask2 <- wind_mask2 / wind_mask2
plot(wind_mask2, main = 'Wind mask alt method', legend = FALSE)

compareRaster(wind_mask_rast, wind_mask2, values = TRUE)

## [1] TRUE

The two masks are identical. Let’s compare to what I got in ArcMap:

wind_key <- raster('hw3_arc_outputs/wind_mask.tif')
check_rast(wind_mask_rast, wind_key)

##   good diff_tot celldiff cellcount
## 1 TRUE        0        0         0

check_rast(wind_mask2, wind_key)

##   good diff_tot celldiff cellcount
## 1 TRUE        0        0         0

And write one to output:

### use compareRaster to compare aspects of the rasters; all.equal is a wrapper
### around compareRaster that compares all aspects of the rasters.
all.equal(wind_mask_rast, wind_mask2) 

## [1] TRUE

writeRaster(wind_mask_rast, 'output/wind_mask.tif', overwrite = TRUE)

Create road mask

Let’s read in the roads layer, and create a mask from that (remembering we need to project it to a new CRS, using sf::st_transform(), and using our ca_alb CRS we borrowed from the County layer). The raster package doesn’t have a nice Euclidean Distance function, and the distance and distanceToPoints functions like an old-school spatial object (sp package), rather than an sf object. So here we’ll use a different method:

• use st_transform() to project the roads into CA Teale Albers.
• feed that into st_buffer and set a distance of 7,500 m.
• collapse all the results into a single feature (we don’t care about the road details, just whether there is a road or not)
• Note, this is a little slow. Let’s check to see if we already have a mask, and if so, skip rerunning it.
• Let’s save it as a geopackage rather than a shapefile!

(We could do this same method in ArcMap too, using different tools instead of Euclidean Distance and Reclassify).

road_mask_file <- 'output/road_mask.gpkg'

if(!file.exists(road_mask_file)) {
  # st_layers('data/HW3.gdb')
  roads_raw_sf <- read_sf(dsn = 'data/HW3.gdb', layer = 'CenterLines2008')
  # st_crs(roads_raw_sf) ### note +units=us-ft instead of meters
  roads_project_sf <- roads_raw_sf %>%
    st_transform(ca_alb)
  
  ### Now, buffer it, union the results (flatten into a single
  ### feature), then intersect it with the County (not necessary,
  ### but why not, trims off the buffer that extends beyond the county line)
  roads_mask_sf <- roads_project_sf %>%
    st_buffer(dist = 7500) %>%
    st_union(by_feature = FALSE) %>%
    st_intersection(county_sf)
  
  write_sf(roads_mask_sf, road_mask_file, delete_layer = TRUE)
} else {
  roads_mask_sf <- read_sf(road_mask_file)
}

plot(roads_mask_sf, main = 'Roads mask', legend = FALSE)

But we can also rasterize this. When putting all the masks together, having them all as rasters can give us some nice advantages. So let’s just do that here while we’re at it. We’ll use the fasterize::fasterize() function to quickly take the sf object and turn it to a raster.

Why fasterize() and not the raster::rasterize()? a) it’s faster, by a long shot. b) rasterize() sometimes chokes when rasterizing polygons that have holes in them, leaving weird artifacts in the resulting raster that will screw up your analysis.

roads_mask_rast <- fasterize::fasterize(roads_mask_sf, base_rast)

plot(roads_mask_rast, 
     main = 'Roads mask as raster', legend = FALSE)

### The vector roads mask has different extensions, so we can save the raster
### as a tif without risk of overwriting the vector files.
writeRaster(roads_mask_rast, 'output/roads_mask.tif', overwrite = TRUE)

Let’s compare to what I got in ArcMap:

roads_key <- raster('hw3_arc_outputs/road_mask.tif')
check_rast(roads_mask_rast, roads_key)

##   good diff_tot    celldiff cellcount
## 1 TRUE      234 0.001627781         0

plot_diffs(roads_mask_rast, roads_key, title = 'Comparing road rast')

### Create urban mask

We can do this process in a similar way we did the roads, using st_buffer to create the exclusion area (rather than the inclusion area), and then subtracting that from the County polygon.

• Union the two Urban layers.
• Feed the urban polygons into st_buffer and set a distance of 1609 m (~1 mile)
• Subtract that from the County polygon using st_difference().
• Let’s save it as a geopackage rather than a shapefile.

Note that I encountered an error:

Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)) : 
   Evaluation error: TopologyException: Input geom 1 is invalid: 
   Ring Self-intersection at or near point ...

These are caused by tiny mistakes in the shapefile that are pretty much impossible to spot visually, where two polygon vertices land right on top of each other. There are some tests built into sf: st_is_valid(), st_is_simple() that will return TRUE or FALSE for each feature in the polygon set. We can use these to find the errors… but sometimes our manipulations (st_union(), st_intersect() etc) cause new glitches.

An easy (but not intuitive) trick can help fix invalid geometries: adding a zero-distance buffer (st_buffer(x, dist = 0)).

urban_mask_file <- 'output/urban_mask.gpkg'

if(!file.exists(urban_mask_file)) {
  # st_layers('data/HW3.gdb')
  urban_2001_sf <- read_sf(dsn = 'data/HW3.gdb', layer = 'Urban2001')
  urban_2020_sf <- read_sf(dsn = 'data/HW3.gdb', layer = 'Urban2020')
  
  # st_is_valid(urban_2001_sf) ### these all look good
  # st_is_valid(urban_2020_sf) ### looks like some invalid geometries in there
  ### can we get around the error by putting a zero buffer on the urban_2020 layer?
  urban_2020_fix <- urban_2020_sf %>% st_buffer(dist = 0)
  # st_is_valid(urban_2020_fix) ### all are valid now!

  ### union the two layers, then make a buffer, then flatten that
  ### using union with only one argument
  urban_buffer_sf <- st_union(urban_2001_sf, urban_2020_fix) %>%
    st_buffer(dist = 1609) %>%
    st_union()
  
  ### now subtract that from the County layer.
  urban_mask_sf <- st_difference(county_sf, urban_buffer_sf)
  
  write_sf(urban_mask_sf, urban_mask_file, delete_layer = TRUE)
} else {
  urban_mask_sf <- read_sf(urban_mask_file)
}

plot(urban_mask_sf %>% select('NAME'), 
     main = 'Urban mask', legend = FALSE)

Again, let’s also make a raster version for later.

urban_mask_rast <- fasterize::fasterize(urban_mask_sf, base_rast)

plot(urban_mask_rast, main = 'Urban mask as raster', legend = FALSE)

writeRaster(urban_mask_rast, 'output/urban_mask.tif', 
            overwrite = TRUE)

Let’s compare to what I got in ArcMap:

urban_key <- raster('hw3_arc_outputs/urban_mask.tif')
check_rast(urban_mask_rast, urban_key)

##    good diff_tot   celldiff cellcount
## 1 FALSE     1985 0.01443668         0

plot_diffs(urban_mask_rast, urban_key, 
           title = 'Comparing urban mask buffer method')

# st_layers('data/HW3.gdb')
urban_2001_sf <- read_sf(dsn = 'data/HW3.gdb', layer = 'Urban2001')
urban_2020_sf <- read_sf(dsn = 'data/HW3.gdb', layer = 'Urban2020')
urban_2020_fix <- urban_2020_sf %>% st_buffer(dist = 0)

### union the two layers to a single feature.  The
### fasterize() function requires it to be a MULTIPOLYGON geometry,
### but for some reason st_union creates a generic geometry, so
### I'll use st_cast() to force it back to MULTIPOLYGON
urban_sf <- st_union(urban_2001_sf, urban_2020_fix) %>%
  st_cast('MULTIPOLYGON')

urban_rast <- fasterize(urban_sf, base_rast)

### Now we can use the raster::distance() function to find the
### distance away from the urban areas
urban_dist_rast <- raster::distance(urban_rast)

plot(urban_dist_rast, main = 'Distance from urban areas', legend = FALSE)

### Now, we can reclassify it (though I'll use the indexing version
### because it just seems cleaner to me).  First create a copy
urban_mask2 <- urban_dist_rast
values(urban_mask2)[values(urban_mask2) <= 1609] <- NA
urban_mask2 <- urban_mask2 / urban_mask2

### Finally, let's mask it with the county 
urban_mask2 <- urban_mask2 %>%
  mask(county_sf)
writeRaster(urban_mask2, 'output/urban_mask2.tif', 
            overwrite = TRUE)

plot(urban_mask2, 
     main = 'Urban mask Euc Dist method', legend = FALSE)

all.equal(urban_mask2, urban_mask_rast) ### FALSE!  not equal.

## [1] FALSE

# compareRaster(urban_mask2, urban_mask_rast, values = TRUE) ### error: not all values equal
plot(urban_mask2, col = 'red', legend = FALSE)
plot(urban_mask_rast, col = 'blue', legend = FALSE, add = TRUE)

urban_key <- raster('hw3_arc_outputs/urban_mask.tif')
check_rast(urban_mask2, urban_key)

##    good diff_tot    celldiff cellcount
## 1 FALSE      515 0.003745536         0

plot_diffs(urban_mask2, urban_key, 
           title = 'Comparing urban mask Euc Dist method')

Create public mask

Again, this builds on what we’ve done with previous layers. We can filter the parcels to only include the public parcels, then subtract those from the county. As with the roads, we’ll need to reproject the layer (at some point) to the CRS of California Teale Albers.

public_mask_file <- 'output/public_mask.gpkg'

if(!file.exists(public_mask_file)) {
  # st_layers('data/HW3.gdb')
  parcels_sf <- read_sf(dsn = 'data/HW3.gdb', layer = 'Parcels')

  ### transform the CRS, then filter for the public use codes.
  ### The same text ordering trick Frew used with the Select tool
  ### works in R (and any other language as far as I know) - so
  ### '8000' is less than '90'.  Let's also flatten it with st_union
  ### so it is easier for the st_difference() function.
  parcels_public <- parcels_sf %>%
    st_transform(ca_alb) %>%
    filter(USECODE < '9000' & USECODE >= '6000') %>%
    st_union(by_feature = FALSE)
           
  ### Subtract that from the County layer.
  public_mask_sf <- st_difference(county_sf, parcels_public)
  
  write_sf(public_mask_sf, public_mask_file, delete_layer = TRUE)
} else {
  public_mask_sf <- read_sf(public_mask_file)
}

plot(public_mask_sf %>% select('NAME'), 
     main = 'Public mask', legend = FALSE)

And yet again, let’s also make a raster version for later.

public_mask_rast <- fasterize::fasterize(public_mask_sf, base_rast)

plot(public_mask_rast, main = 'Public mask as raster', legend = FALSE)

writeRaster(public_mask_rast, 'output/public_mask.tif', 
            overwrite = TRUE)

Create fire mask

This is nearly identical to the Public mask. On my first attempt, I noticed an error that causes things to choke:

Evaluation error: TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point 50040.80700000003 -400081.27190000005 at 50040.80700000003 -400081.27190000005.

To get around this, I dissolve the interior boundaries (using st_union as a “dissolve” tool) of the two fire layers before joining them. Sometimes complicated polygon features have conflicts when you join them or subtract them; ArcMap is good at working around these conflicts, but R is not so good. But by dissolving a bunch of polygons into one, we can sometimes avoid creating those conflicts. If you need to keep track of the attributes, this workaround will not work for you.

fire_mask_file <- 'output/fire_mask.gpkg'

if(!file.exists(fire_mask_file)) {
  # st_layers('data/HW3.gdb')
  fire_lra_sf <- read_sf(dsn = 'data/HW3.gdb', layer = 'FireLRA')
  fire_sra_sf <- read_sf(dsn = 'data/HW3.gdb', layer = 'FireSRA')

  ### Filter the SRA to the proper hazard codes, and then union
  ### the two layers.  Again, let's also flatten it with st_union
  ### so it is easier for the st_difference() function.
  fire_sra_filtered_sf <- fire_sra_sf %>%
    filter(HAZ_CODE == 3) %>%
    st_union(by_feature = FALSE) ### flatten it
  
  fire_joined_sf <- fire_lra_sf %>%
    st_union(by_feature = FALSE) %>% ### flatten it
    st_union(fire_sra_filtered_sf)   ### 
           
  ### Subtract that from the County layer.
  fire_mask_sf <- st_difference(county_sf, fire_joined_sf)
  
  write_sf(fire_mask_sf, fire_mask_file, delete_layer = TRUE)
} else {
  fire_mask_sf <- read_sf(fire_mask_file)
}

plot(fire_mask_sf %>% select('NAME'), 
     main = 'Fire mask', legend = FALSE)

And the raster version.

fire_mask_rast <- fasterize::fasterize(fire_mask_sf, base_rast)

plot(fire_mask_rast, main = 'Fire mask as raster', legend = FALSE)

writeRaster(fire_mask_rast, 'output/fire_mask.tif', 
            overwrite = TRUE)

Create airport mask

This is nearly identical to the Urban mask - I’ll do this with a buffer again.

airport_mask_file <- 'output/airport_mask.gpkg'

if(!file.exists(airport_mask_file)) {
  # st_layers('data/HW3.gdb')
  airport_sf <- read_sf(dsn = 'data/HW3.gdb', layer = 'Airports')

  airport_buffer_sf <- airport_sf %>%
    st_buffer(dist = 7500) %>%
    st_union(by_feature = FALSE)
  
  ### Subtract that from the County layer.
  airport_mask_sf <- st_difference(county_sf, airport_buffer_sf)
  
  write_sf(airport_mask_sf, airport_mask_file, delete_layer = TRUE)
} else {
  airport_mask_sf <- read_sf(airport_mask_file)
}

plot(airport_mask_sf %>% select('NAME'), 
     main = 'airport mask', legend = FALSE)

And the raster version.

airport_mask_rast <- fasterize::fasterize(airport_mask_sf, base_rast)

plot(airport_mask_rast, main = 'Airport mask as raster', legend = FALSE)

writeRaster(airport_mask_rast, 'output/airport_mask.tif', 
            overwrite = TRUE)

airport_key <- raster('hw3_arc_outputs/airport_mask.tif')
check_rast(airport_mask_rast, airport_key)

##    good diff_tot    celldiff cellcount
## 1 FALSE      546 0.002355438         0

plot_diffs(airport_mask_rast, airport_key, 
           title = 'Comparing airport mask')